#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <cstdlib>
#include "txtpara.h"
#include "KDindex.h"
#include "directsampling.h"
#include "tplate.h"
#include <boost/random.hpp>
#include <ctime>
using namespace std;
double inline normalize(double y, double min, double max, double maxdiff)
{
    return (y-min)/(maxdiff)*255;
}
// Generating normally distributed random numbers
double sampleNormal(double mean, double sigma)
{
    //Establish Mersenne twister random number generator
    //Set seed using UNIX time
    static boost::mt19937 rng(static_cast < unsigned >(std::time(0)));
    //Select Gaussian random distribution
    boost::normal_distribution < double >norm_dist(mean, sigma);
    //Use the form of function to generate a random number generator
    boost::variate_generator < boost::mt19937 &,
        boost::normal_distribution < double >>normal_sampler(rng,norm_dist);
    return normal_sampler();
}

int main(int argc, char *argv[])
{
    if (argc != 5)
    {
        cout << "Usage:\n" << argv[0] <<
            " param_file training_image_file template_file paramfile"
            << endl;
        return 0;
    }
    /* direct sampling parameters */
    int NR_neighbor;            //Maximum number of neighbors
    int NR_try;                    //Maximum number of search
    float tolerance;            //tolerance of difference
    int xsize, ysize, zsize;
    double percentage;
    int unitsize = 1;
    int maxdv = 20;
    char *config = argv[1];
    NR_neighbor = getparaINT("NR_neighbor", config);
    tolerance = getparaFLOAT("tolerance", config);
    xsize = getparaINT("xsize", config);
    ysize = getparaINT("ysize", config);
    zsize = getparaINT("zsize", config);
    maxdv = getparaINT("maxdv", config);
    percentage = getparaFLOAT("percent", config);
    std::vector < coord > sample;
    //Step 1.  Generate sample points for the first time
    //>>>>>1.1 calculating Number of points
    int tol_sample = xsize * ysize * zsize * percentage;
    int N_sample = xsize * ysize * zsize * percentage;
    int interval = round(sqrt(xsize * ysize / N_sample));
    if ((xsize/(interval))*(ysize/(interval))>N_sample)
        interval++;
    //>>>>>1.2 generate random path
    srand(clock());
    std::vector < int >ids;
    int zpos = 0;
    int s_no = 0;
    for (int xpos = 0; xpos < xsize; xpos += interval)
    {
        for (int ypos = 0; ypos < ysize; ypos += interval)
        {
            coord s;
            s.x = xpos;
            s.y = ypos;
            s.z = zpos;
            sample.push_back(s);
            ids.push_back(s_no);
            s_no++;
        }
    }
    N_sample = sample.size();
    for(int xpos=0;xpos<xsize;xpos++)
    {
        for(int ypos=0;ypos<ysize;ypos++){
            if(((xpos % interval)==0) && ((ypos % interval)==0))
                continue;
            coord s;
            s.x = xpos;
            s.y = ypos;
            s.z = zpos;
            sample.push_back(s);
            ids.push_back(s_no);
            s_no++;
        }
    }
    for(int num=0;num<5;num++){
        for(int i=N_sample;i<tol_sample;i++)
        {
            long rndnum=0;
            for(int j=0;j<(xsize*ysize/RAND_MAX+10);j++){
                rndnum+=rand();
            }
            int ranid = N_sample+rndnum % (xsize*ysize-N_sample);
            int tmp = ids[i];
            ids[i] = ids[ranid];
            ids[ranid] = tmp;
            //Take the rest of the random samples from the rest
        }
    }
    //Random order
    for (int s_no = 0; s_no < tol_sample; s_no++)
    {
        int ranid = rand() % tol_sample;
        int tmp = ids[s_no];
        ids[s_no] = ids[ranid];
        ids[ranid] = tmp;
    }
    N_sample=tol_sample;
    //Step 2. Perform the simulation for first fixed percentage of cells
    Index < L2 < float >>*R = createIndex();
    float *ti;
    int tilength, tiwidth, tiheight;
    int start_id = 0;
    int level = 0;
    float *target = new float[int (xsize * ysize * zsize)];
    float *target_ = new float[int (xsize * ysize * zsize)];
    for (int xpos = 0; xpos < xsize; xpos++)
    {
        for (int ypos = 0; ypos < ysize; ypos++)
        {
            target[xpos * ysize + ypos] = -1;
        }
    }
    for (int xpos = 0; xpos < xsize; xpos++)
    {
        for (int ypos = 0; ypos < ysize; ypos++)
        {
            target_[xpos * ysize + ypos] = -1;
        }
    }
    load_ti_3D(&ti, argv[2], &tilength, &tiwidth, &tiheight);
    NR_try =
        int (getparaFLOAT("Fraction", config) * tilength * tiwidth *
                tiheight);
    Simulation_cont(R, ti, NR_neighbor, NR_try, tolerance, sample, ids,
            target, N_sample, tilength, tiwidth, tiheight,
            unitsize, start_id, level, xsize, ysize, zsize);
    //Step 3 Interpolation using regression
    //1. Read in template file tp.par
    tplate tp;
    load_template(&tp, argv[3]);
    data_event_cont *dv = new data_event_cont;
    dv->tpos = new int[tp.tpsize];
    dv->pos = new int[2 * tp.tpsize];
    dv->value = new float[tp.tpsize];

    //2. Read in model file model.txt
    /*
     * load theta for estimation
     * */
    ifstream param;
    param.open(argv[4]);
    int nr_param;
    int add_param=3; //The maximum, minimum, and maximum minus the minimum
    int tpsize;
    param >> nr_param;
    param >> tpsize;
    //int effect_size = tpsize;
    float *theta = new float[tpsize * (nr_param+add_param)];
    for (int i = 0; i < tpsize; i++)
    {
        for (int j = 0; j < nr_param; j++)
        {
            param >> theta[i * (nr_param+add_param) + j];
        }
        theta[i*(nr_param+add_param)+nr_param]=theta[i*(nr_param+add_param)];
        theta[i*(nr_param+add_param)+nr_param+1]=theta[i*(nr_param+add_param)]+255.0*theta[i*(nr_param+add_param)+1];
        theta[i*(nr_param+add_param)+nr_param+2]=theta[i*(nr_param+add_param)+nr_param+1]-theta[i*(nr_param+add_param)+nr_param];
    }
    param.close();
    //3. Estimate non simulated points
    for (int xpos = 0; xpos < xsize; xpos++)
    {
        for (int ypos = 0; ypos < ysize; ypos++)
        {
            if ((target[xpos * ysize + ypos] - (-1)) < 1e-6)
            {
                //Not simulated
                //Find the neighbors
                get_dv_cont(dv, &tp, target, xsize, ysize, xpos, ypos,
                        maxdv);
                if (dv->dvsize > 0)    //It has neighbors
                {
                    //Regression results were calculated using each neighbor
                    //Average regression results                   
                    float y;
                    float y_sum = 0;
                    float y_avg;
                    float q;
                    float var_sum = 0;
                    float var_avg;
                    int k = 0;
                    /*
                     * Solve the sum of y values y_sum and the sum of variance var_sum in all conditional data y=ax+b
                     * */
                    for (int i = 0; i < dv->dvsize; i++)
                    {
                        if (0.01 > theta[(dv->tpos[i]) * (nr_param+add_param) + 3])
                            //The pvalue is too large and has no statistical significance. The later ones do not need to be read in
                        {
                            y = theta[(dv->tpos[i]) * (nr_param+add_param)] + theta[(dv->tpos[i]) * (nr_param+add_param) + 1] * (dv->value[i]);    //Obtain the value of each regression y of conditional data
                            y = normalize(y,theta[(dv->tpos[i])*(nr_param+add_param)+nr_param],theta[(dv->tpos[i])*(nr_param+add_param)+nr_param+1],theta[(dv->tpos[i])*(nr_param+add_param)+nr_param+2]);
                            //The number obtained here is equivalent to compressing the result, so you need to decompress the result to 0-255
                            q = theta[(dv->tpos[i]) * (nr_param+add_param) + 2];//Assign the square of the standard deviation to q                                                      
                            q = normalize(q,theta[(dv->tpos[i])*(nr_param+add_param)+nr_param],theta[(dv->tpos[i])*(nr_param+add_param)+nr_param+1],theta[(dv->tpos[i])*(nr_param+add_param)+nr_param+2]);
                            q=q*q;
                            y_sum += y;    //Find the y sum of all conditional data
                            var_sum += q;    //Sum the variance of all conditional data
                            k++;    //Find the number of conditional data                              
                        }
                    }
                    y_avg = y_sum / k;    //Obtain the mean value of conditional data y_avg
                    var_avg = sqrt(var_sum) / k;    //Get sigma                         
                    {
                        target_[xpos * ysize + ypos] = sampleNormal(y_avg, var_avg);    //target_ temporary storage simulation point
                        if (target_[xpos * ysize + ypos] < y_avg-var_avg)
                            target_[xpos*ysize+ypos] = y_avg-var_avg;
                        if (target_[xpos*ysize+ypos] < 0) 
                            target_[xpos*ysize+ypos] = 0;
                        if (target_[xpos * ysize + ypos] > y_avg+var_avg)
                            target_[xpos*ysize+ypos] = y_avg+var_avg;
                        if (target_[xpos*ysize+ypos] > 255) 
                            target_[xpos*ysize+ypos] = 255;
                    }
                }
                else
                {
                    //The result is a random value between the maximum and minimum
                    target_[xpos * ysize + ypos] = (float) (rand() % 255);
                }
            }
        }
    }

    for (int xpos = 0; xpos < xsize; xpos++)
    {
        for (int ypos = 0; ypos < ysize; ypos++)
        {

            if ((target[xpos * ysize + ypos] - (-1)) < 1e-6)
            {
                target[xpos * ysize + ypos] = target_[xpos * ysize + ypos];
            }
        }
    }
    //Step 4 Output the simulation results
    for (int xpos = 0; xpos < xsize; xpos++)
    {
        for (int ypos = 0; ypos < ysize; ypos++)
        {
            cout << target[xpos * ysize + ypos] << " ";
        }
        cout << endl;
    }
    return 0;
}
